Sampling and homology via bottlenecks¶

The paper Sampling and homology via bottlenecks by Di Rocco et. al. introduces a novel algorithm for producing a provably dense sampling of a smooth and compact algebraic variety where the density is determined by the smallest bottleneck. Using the sample, the zeroth and first homology of the variety can be recovered using a modified Vietoris-Rips construction. This notebook implements the algorithm for sampling and homology computation for the case of complete intersections and illustrates it using a curve in $\Bbb R^2$.

In [17]:
using HomotopyContinuation, DynamicPolynomials, LinearAlgebra, IterTools, Random

Random.seed!(1)

n=5 # ambient dimension
@polyvar x[1:n] y[1:n] p[1:n] gamma[1:n]

singular_cubic = x[1]*(x[2]^2+x[3]^2+x[4]^2+x[5]^2)-(x[2]^3+x[3]^3+x[4]^3-(1/2)*x[5]^3);
smoothing = singular_cubic-(1/4)*x[1]^3;

deg_three_generic = rand((-1,1)) + rand((-1,1))*x[1] + rand((-1,1))*x[1]^2 + rand((-1,1))*x[1]^3 +  rand((-1,1))*x[2] + rand((-1,1))*x[1]*x[2] + rand((-1,1))*x[1]^2*x[2] +  rand((-1,1))*x[2]^2 + rand((-1,1))*x[1]*x[2]^2 + rand((-1,1))*x[2]^3 +  rand((-1,1))*x[3] + rand((-1,1))*x[1]*x[3] + rand((-1,1))*x[1]^2*x[3] +  rand((-1,1))*x[2]*x[3] + rand((-1,1))*x[1]*x[2]*x[3] + rand((-1,1))*x[2]^2*x[3] +  rand((-1,1))*x[3]^2 + rand((-1,1))*x[1]*x[3]^2 + rand((-1,1))*x[2]*x[3]^2 +  rand((-1,1))*x[3]^3 + rand((-1,1))*x[4] + rand((-1,1))*x[1]*x[4] +  rand((-1,1))*x[1]^2*x[4] + rand((-1,1))*x[2]*x[4] + rand((-1,1))*x[1]*x[2]*x[4] +  rand((-1,1))*x[2]^2*x[4] + rand((-1,1))*x[3]*x[4] + rand((-1,1))*x[1]*x[3]*x[4] +  rand((-1,1))*x[2]*x[3]*x[4] + rand((-1,1))*x[3]^2*x[4] + rand((-1,1))*x[4]^2 +  rand((-1,1))*x[1]*x[4]^2 + rand((-1,1))*x[2]*x[4]^2 + rand((-1,1))*x[3]*x[4]^2 +  rand((-1,1))*x[4]^3 + rand((-1,1))*x[5] + rand((-1,1))*x[1]*x[5] +  rand((-1,1))*x[1]^2*x[5] + rand((-1,1))*x[2]*x[5] + rand((-1,1))*x[1]*x[2]*x[5] +  rand((-1,1))*x[2]^2*x[5] + rand((-1,1))*x[3]*x[5] + rand((-1,1))*x[1]*x[3]*x[5] +  rand((-1,1))*x[2]*x[3]*x[5] + rand((-1,1))*x[3]^2*x[5] + rand((-1,1))*x[4]*x[5] +  rand((-1,1))*x[1]*x[4]*x[5] + rand((-1,1))*x[2]*x[4]*x[5] +  rand((-1,1))*x[3]*x[4]*x[5] + rand((-1,1))*x[4]^2*x[5] + rand((-1,1))*x[5]^2 +  rand((-1,1))*x[1]*x[5]^2 + rand((-1,1))*x[2]*x[5]^2 + rand((-1,1))*x[3]*x[5]^2 +  rand((-1,1))*x[4]*x[5]^2 + rand((-1,1))*x[5]^3;
deg_two_generic = rand((-1,1)) + rand((-1,1))*x[1] + rand((-1,1))*x[1]^2 + rand((-1,1))*x[2] +  rand((-1,1))*x[1]*x[2] + rand((-1,1))*x[2]^2 + rand((-1,1))*x[3] +  rand((-1,1))*x[1]*x[3] + rand((-1,1))*x[2]*x[3] + rand((-1,1))*x[3]^2 +  rand((-1,1))*x[4] + rand((-1,1))*x[1]*x[4] + rand((-1,1))*x[2]*x[4] +  rand((-1,1))*x[3]*x[4] + rand((-1,1))*x[4]^2 + rand((-1,1))*x[5] +  rand((-1,1))*x[1]*x[5] + rand((-1,1))*x[2]*x[5] + rand((-1,1))*x[3]*x[5] +  rand((-1,1))*x[4]*x[5] + rand((-1,1))*x[5]^2;

epsilon = 1/100

F = [
    smoothing+epsilon*deg_three_generic,
    x[1]^2+x[2]^2+x[3]^2+x[4]^2+x[5]^2-1+epsilon*deg_two_generic
]

d=length(F) # codimension of variety
k = n-d # dimension of variety
@polyvar lambda[1:d] mu[1:d]; # lagrange multipliers
In [18]:
# Compute the bottlenecks

grad = differentiate(F, x)
G = subs(F, x => y)
grady = subs(grad, x => y)

system = [F; G; map(j -> x[j]-y[j]-dot(lambda, grad[:, j]), 1:n); map(j -> x[j]-y[j]-dot(mu, grady[:, j]), 1:n)]
result = solve(system, start_system = :total_degree)

Tracking 2125764 paths... 100%|████████████████████████████████████████| Time: 1:53:40
  # paths tracked:                  2125764
  # non-singular solutions (real):  20724 (712)
  # singular solutions (real):      5721 (12)
  # total solutions (real):         26445 (724)
Out[18]:
Result{Vector{ComplexF64}} with 23856 solutions
===============================================
• 20724 non-singular solutions (712 real)
• 3132 singular solutions (12 real)
• 2125764 paths tracked
• random seed: 835247
• multiplicity table of singular solutions:
┌───────┬───────┬────────┬────────────┐
│ mult. │ total │ # real │ # non-real │
├───────┼───────┼────────┼────────────┤
│   1   │  543  │   12   │    531     │
│   2   │ 2589  │   0    │    2589    │
└───────┴───────┴────────┴────────────┘
In [19]:
# Choose the size of the grid as the smallest bottleneck

bottlenecks = map(s -> (s[1:n], s[n+1:2*n]), real_solutions(nonsingular(result)))
bn_lengths = sort!(map(b -> (norm(b[1]-b[2]), b), bottlenecks), by = a -> a[1])
δ = bn_lengths[1][1]/2 # grid size
Out[19]:
0.09505817482972509
In [20]:
# Compute the bounding box by computing the EDD starting from the center of the largest bottleneck

q = bn_lengths[end][2][1] + (bn_lengths[end][2][2]-bn_lengths[end][2][1])/2 + randn(Float64, n)*10^(-10)
system = [F; map(j -> x[j]-q[j]-dot(lambda, grad[:, j]), 1:n)]
result = solve(system, start_system = :polyhedral)
Out[20]:
Result{Vector{ComplexF64}} with 156 solutions
=============================================
• 156 non-singular solutions (24 real)
• 0 singular solutions (0 real)
• 156 paths tracked
• random seed: 104798
In [21]:
# Extract farthest point from q to X and use as box length

critical_points = sort!(map(c -> (norm(c[1:n]-q), c[1:n]), real_solutions(nonsingular(result))), by = a -> a[1])
b = critical_points[end][1]
indices = [i for i in -b:δ:b];
In [22]:
# Compute basic sample

samples = []
counter = 0

start_time = time_ns()
for s in IterTools.subsets(1:n, k)
    Ft = [F; map(i -> x[s[i]]-p[i]-q[s[i]], 1:k)]
    p₀ = randn(ComplexF64, k)
    F_p₀ = subs(Ft, p[1:k] => p₀)
    result_p₀ = solve(F_p₀)
    S_p₀ = solutions(result_p₀)
    
    # Construct the PathTracker
    tracker = HomotopyContinuation.pathtracker(Ft; parameters=p[1:k], generic_parameters=p₀)
    for p1 in Iterators.product(map(j-> 1:length(indices), s)...)
        counter += length(S_p₀)
        for s1 in S_p₀
            result = track(tracker, s1; target_parameters=map(j -> indices[p1[j]], 1:k))
            # check that the tracking was successfull
            if is_success(result) && is_real(result)
                push!(samples, real(solution(result)))
            end
        end
    end
end
time_basic_sample = time_ns()-start_time;
In [23]:
# Compute extra sample

extra_samples = []
extra_counter = 0

start_time = time_ns()
for l in 1:k-1
    for s in IterTools.subsets(1:n, l)
        Ft = [F; map(i -> x[s[i]]-p[i]-q[s[i]], 1:l)] 
        grad = differentiate(Ft, x)
        system = [Ft; map(j -> x[j]-y[j]-dot(gamma[1:n-k+l], grad[:, j]), 1:n)]
        
        p₀ = randn(ComplexF64, n+l)
        F_p₀ = subs(system, [y; p[1:l]] => p₀)
        result_p₀ = solve(F_p₀)
        S_p₀ = solutions(result_p₀)

        # Construct the PathTracker
        tracker = HomotopyContinuation.pathtracker(system; parameters=[y; p[1:l]], generic_parameters=p₀)
        for p1 in Iterators.product(map(j-> 1:length(indices), s)...)
            extra_counter += length(S_p₀)
            for s1 in S_p₀
                result = track(tracker, s1; target_parameters=[randn(Float64, n); map(j -> indices[p1[j]], 1:l)])
                # check that the tracking was successfull
                if is_success(result) && is_real(result)
                    push!(extra_samples, real(solution(result))[1:n])
                end
            end
        end
    end
end
time_extra_sample = time_ns()-start_time;

Tracking 1458 paths... 100%|████████████████████████████████████████| Time: 0:00:00
  # paths tracked:                  1458
  # non-singular solutions (real):  24 (0)
  # singular solutions (real):      0 (0)
  # total solutions (real):         24 (0)
In [24]:
# Summary

println("grid size: ", δ)
println("length of bounding box: ", b)
println("number of tracked paths: ", counter+extra_counter)
println("total time: ", (time_basic_sample+time_extra_sample)*10^(-9), " s")
println("|E_δ|: " , length(samples), ", |E'_δ|: ", length(extra_samples))
println("total number of samples: ", length(samples)+length(extra_samples))

grid size: 0.09505817482972509
length of bounding box: 1.0137680752023333
number of tracked paths: 762300
total time: 83.85367645500004 s
|E_δ|: 69154, |E'_δ|: 10841
total number of samples: 79995

Export samples¶

In [25]:
using DelimitedFiles
writedlm( "pointcloud.csv",  samples, ',')

Plotting the samples¶

In [26]:
using Plots
plotlyjs()

ArgumentError: Package PlotlyJS not found in current path:
- Run `import Pkg; Pkg.add("PlotlyJS")` to install the PlotlyJS package.


Stacktrace:
  [1] require(into::Module, mod::Symbol)
    @ Base ./loading.jl:967
  [2] top-level scope
    @ ~/.julia/packages/Plots/hyS17/src/backends.jl:491
  [3] eval
    @ ./boot.jl:373 [inlined]
  [4] _initialize_backend(pkg::Plots.PlotlyJSBackend)
    @ Plots ~/.julia/packages/Plots/hyS17/src/backends.jl:490
  [5] backend(pkg::Plots.PlotlyJSBackend)
    @ Plots ~/.julia/packages/Plots/hyS17/src/backends.jl:174
  [6] #plotlyjs#242
    @ ~/.julia/packages/Plots/hyS17/src/backends.jl:31 [inlined]
  [7] plotlyjs()
    @ Plots ~/.julia/packages/Plots/hyS17/src/backends.jl:31
  [8] top-level scope
    @ In[26]:2
  [9] eval
    @ ./boot.jl:373 [inlined]
 [10] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
    @ Base ./loading.jl:1196
In [27]:
# Plot the sampled points

S = reduce(hcat, vcat(samples, extra_samples))
if n < 3
    scatter(S[1,:], S[2,:])
else
    scatter(S[1,:], S[2,:], S[3,:])
end
Out[27]:

Homology¶

In [28]:
using Eirene
In [29]:
#Compute distance matrix

ϵ = 2*δ
D = zeros((length(samples), length(samples))) #Initialize distance matrix

neighbour_lists = []
for i in 1:length(samples)
    push!(neighbour_lists, [])
end
candidate_lists = []

for i in 1:length(samples)
    candidate_list = []
    for j in (i+1):length(samples)
        dist = norm(samples[i]-samples[j])
        if dist < sqrt(8)*δ
            if dist < ϵ
                push!(neighbour_lists[i], j)
                push!(neighbour_lists[j], i)
            else
                push!(candidate_list, j)
            end
        end
        D[i, j] = dist
        D[j, i] = dist
    end
    push!(candidate_lists, candidate_list)
end

OutOfMemoryError()

Stacktrace:
 [1] Array
   @ ./boot.jl:459 [inlined]
 [2] Array
   @ ./boot.jl:467 [inlined]
 [3] zeros
   @ ./array.jl:525 [inlined]
 [4] zeros(dims::Tuple{Int64, Int64})
   @ Base ./array.jl:522
 [5] top-level scope
   @ In[29]:4
 [6] eval
   @ ./boot.jl:373 [inlined]
 [7] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1196
In [30]:
# Modify distance matrix to add edges in VR-complex

DM = deepcopy(D) #Modified distance matrix

thresh = 2*δ - 10^(-10)
for i in 1:length(samples)
    for j in candidate_lists[i]
        for k in neighbour_lists[j]
            if D[i, k] < ϵ && D[j, k] < ϵ
                DM[i, j] = thresh
                DM[j, i] = thresh
                break
            end
        end
    end
end

UndefVarError: D not defined

Stacktrace:
 [1] top-level scope
   @ In[30]:3
 [2] eval
   @ ./boot.jl:373 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1196
In [31]:
# Homology computation using Eirene using VR-complex with distance ϵ
C = eirene(DM, model="vr", maxdim=1, minrad=ϵ, maxrad=ϵ);

UndefVarError: DM not defined

Stacktrace:
 [1] top-level scope
   @ In[31]:2
 [2] eval
   @ ./boot.jl:373 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1196
In [32]:
println("0-th Betti number: ", betticurve(C, dim=0)[1, 2])
println("1-th Betti number: ", betticurve(C, dim=1)[1, 2])

UndefVarError: C not defined

Stacktrace:
 [1] top-level scope
   @ In[32]:1
 [2] eval
   @ ./boot.jl:373 [inlined]
 [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)
   @ Base ./loading.jl:1196